No evidence of neural feature-specific pre-activation during the prediction of an upcoming stimulus

A reanalysis of Demarchi et al. (Nature Communications, 2019)

Authors

Affiliation

Abdoun Oussama

Centre de Rercherches en Neurosciences de Lyon

Todorov Dimitrii

Poublan-couzardot Arnaud

Tirou Coumarane

Lutz Antoine

Vernet Marine

Quentin Romain

README

Background information

See the associated Github project home page.

The present R Quarto document calculates the correlation between decoding accuracy and sequence entropy, perform cluster-based permutation statistical tests and Bayesian statistics on accuracy and correlation time-generalization matrices, and plot all the figures and supplementary figures of the paper.

Instructions for use

1. Obtain the classification data either from Zenodo, or by running the Python scripts on raw MEG data (see the Github project home page for details). Save it in the ./data folder.

2. For a first quick run of the script, we recommend that you set the number of cluster permutations to a low value (e.g. 100) in Prepare / Parameters / Global. Set it to 10,000 for stable results matching the ones reported in the publication.

Prepare

Libraries

library(expm)         # exponentiate matrices

Loading required package: Matrix

Attaching package: 'expm'

The following object is masked from 'package:Matrix':

    expm

library(patchwork)    # combine plots
library(tictoc)       # time executions
library(correlation)
library(abind)        # combine multi-dimensional arrays
library(BayesFactor)

Loading required package: coda

************
Welcome to BayesFactor 0.9.12-4.6. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).

Type BFManual() to open the manual.
************

library(magrittr)     # special pipes
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.0     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.0
✔ purrr     1.0.2     

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand()    masks Matrix::expand()
✖ tidyr::extract()   masks magrittr::extract()
✖ dplyr::filter()    masks stats::filter()
✖ dplyr::lag()       masks stats::lag()
✖ tidyr::pack()      masks Matrix::pack()
✖ purrr::set_names() masks magrittr::set_names()
✖ tidyr::unpack()    masks Matrix::unpack()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(reticulate)
np <- import("numpy")
os <- import("os")
mne <- import("mne")

Parameters

Global

path.root <- "./data/results"
list.subj <- list.dirs(path.root, full.names = F, recursive = F)
n.subj <- length(list.subj)

time_offset <- 670
n_permutations <- 100
save <- FALSE

Plotting

# Temporal resolution in ms (1000 / sampling frequency)
dt <- 10

# Colors
col.pal.base <- rev(pals::brewer.rdbu(10))
col.pal.diff <- rev(pals::brewer.brbg(10))
col.pal.corr <- pals::ocean.curl(100)
col.lines <- "grey40"
col.clusters <- "grey10"

# Breaks & Labels
# --- function
format_labels_acc <- function(labels) {
  labels[seq(2,length(labels),2)] <- ""
  labels <- str_replace(labels, "0.", ".")
}
# --- base
acc.breaks.base <- seq(.21,.29, .005)
acc.labels.base <- format_labels_acc(acc.breaks.base)
# --- diff
acc.breaks.diff <- seq(-.02,.02, .005)
acc.labels.diff <- format_labels_acc(acc.breaks.diff)
# --- correlations
corr.breaks <- seq(-1,1,.2)
corr.labels <- str_replace(corr.breaks, "0.", ".")

# GGPLOT theme
theme_set(theme_minimal())
theme_update(text = element_text(family = "Arial"),
             plot.tag.position = c(0, 0.95),
             plot.tag = element_text(size = 14, hjust = 0, vjust = 1, face = "bold"))

Transition matrices

From Demarchi, Sanchez, and Weisz (2019)

mat.T <- list(RD = matrix(rep(.25, 16), ncol = 4, byrow = T),
              MM = matrix(), MP = matrix(), OR = matrix())

mat.T$MM <- matrix(c(c(.25,0,.37,.38),
                     c(.38,.25,0,.37),
                     c(.37,.38,.25,0),
                     c(0,.37,.38,.25)),
                   ncol = 4, byrow = T)

mat.T$MP <- matrix(c(c(.25,0,.15,.60),
                     c(.60,.25,0,.15),
                     c(.15,.60,.25,0),
                     c(0,.15,.60,.25)),
                   ncol = 4, byrow = T)

mat.T$OR <- matrix(c(c(.25,0,0,.75),
                     c(.75,.25,0,0),
                     c(0,.75,.25,0),
                     c(0,0,.75,.25)),
                   ncol = 4, byrow = T)

covCT <- function(mat.C, mat.T) {
  covar <- 0
  for (i in 1:dim(mat.C)[1]) {
    covar <- covar + mean((mat.C[i,]-mean(mat.C[i,]))*(mat.T[i,]-mean(mat.T[i,])))
  }
  return(covar)
}

Functions

Matrix to df

matrix_to_df2 <- function(array, row.name, col.name, val.name) {
  array %>% 
    as_tibble() %>%
    rownames_to_column(row.name) %>%
    pivot_longer(-row.name, names_to = col.name, values_to = val.name) %>% 
    mutate("{col.name}" := str_remove(!!sym(col.name), "V")) %>% 
    mutate_all(as.numeric) %>% 
    return
}

# !!! USE as.data.frame.table() instead

Batch data loading

load_loop_subj <- function(root = path.root, folder, condition, diag = F, crop=T) {
  
  # Individual scores
  if (folder == "omissions") {
    array.scores <- array(dim = c(141,141,n.subj))
  } else if (folder == "sounds") {
    array.scores <- array(dim = c(33,135,n.subj))
  }
  
  i <- 1
  for (s in list.subj) {
    file <- os$path$join(root, s, folder,
                         paste0("cv_", condition, "_scores.npy"))

    # Proceed to next subject if file doesn't exist
    if (!os$path$isfile(file)) {
      cat(file, " does not exist !!!\n")
      next
    }
    
    # Load data otherwise
    tmp <- np$load(file)
    
    # Average across CV folds if present
    if (length(dim(tmp)) == 3) {
      tmp <- tmp %>% apply(c(2,3), mean) # average across CV folds
    }
      
    array.scores[,,i] <- tmp
    i <- i+1
  }
  
  # Average across participants
  t_train_offset <- case_when(folder == "omissions" ~ time_offset,
                              folder == "sounds" ~ 0)
  t_test_offset <- time_offset
  
  df.mean <- array.scores %>% 
    apply(c(1,2), mean) %>%
    matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "accuracy") %>%
    # --- convert indices to time
    mutate(t_train = dt*(t_train-1) - t_train_offset,
           t_test =  dt*(t_test-1)  - t_test_offset)
  
  # --- crop time if requested (default)
  if (crop) {
    df.mean %<>% filter(t_train>=0, (t_train-333) < dt)
  }
  
  return(list(array.scores, df.mean))
}

Correlations

## Correlation function between two vectors
correlations_vecs <- function(vec1, vec2) {
  cor(as.vector(vec1), as.vector(vec2))
}

# Bind time-generalization matrices from different entropy conditions into a single 3D array
correlations_extract3D <- function(data, idx.lines, idx.subj) {
  arr.tmp <- data[[idx.lines[1],"array"]][[1]][,,idx.subj]
  for (i in 2:length(idx.lines)) {
    arr.tmp <- abind(arr.tmp, data[[idx.lines[i],"array"]][[1]][,,idx.subj], along = 3)
  }
  return(arr.tmp)
}


# Replace NAs by 0 and extreme values by +/- .99
correlations_fixval <- function(arr.corr) {
  arr.corr[is.na(arr.corr)] <- 0
  arr.corr[arr.corr == 1] <- .99
  arr.corr[arr.corr == -1] <- -.99
  return(arr.corr)
}

Cluster permutation stats

get_signif_clust2d <- function(X, n_permutations = 2**12, n_jobs = -1) {
  # run permutaiton tests
  mne$stats$spatio_temporal_cluster_1samp_test(X,
                                               n_permutations = n_permutations,
                                               n_jobs = as.integer(n_jobs),
                                               out_type = "mask",
                                               verbose = FALSE) -> tmp
  names(tmp) <- c("t_obs","clusters","cluster_pv","H0")
  
  # Extract clusters: the final output is a 0/1 matrix, 1 indicating belonging to a significant cluster 
  # --- get ids of significant clusters
  idx.signif <- which(tmp$cluster_pv < .05) 
  # --- aggregate clusters
  signif <- array(0, dim = dim(tmp$t_obs)) 
  for (i in idx.signif) {
    signif <- signif + tmp$clusters[[i]]
  }
  
  return(1*(signif > 0))
}


# Loop through multiple conditions, stored in a dataframe where each row contains a 3D array in an "array" column
# Other columns are appended to the output as descriptors of the condition
cluster_loop_cond <- function(df, H0 = 0, n = 2**12, n_jobs = -1, crop = FALSE) {
  # Initialize significant clusters dataframe
  df.clusters <- tibble()
  # Loop through all conditions
  for (i in 1:nrow(df)) {
    
    # Extract array of subject-level time-generalized accuracy
    X <- df[[i,"array"]][[1]]
    # --- crop to stimulus N-1 if requested
    if (crop) {X <- X[,35:68,]}
    # --- reshape to how MNE expects it
    X <- X %>% np$moveaxis(as.integer(2), as.integer(0)) 
  
    df.clusters %<>% bind_rows(
      get_signif_clust2d(X-H0, n_jobs = n_jobs, n_permutations = n) %>% 
        matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "signif") %>% 
        bind_cols(select(df,-array)[i,])
    )
  }
  
  return(df.clusters)
}

Difference original-reordered

calculate_diff <- function(df, f.map, values_to) {
  df %>% 
    mutate(manip = case_match(manip, "" ~ "original", "_reord" ~ "reordered")) %>% 
    pivot_wider(names_from = manip, values_from = values_to) %>% 
    mutate(diff = f.map(original, reordered, .f=\(x,y){x-y})) %>% 
    pivot_longer(cols = c(original, reordered, diff), names_to = "manip", values_to = values_to)
}

Base plot

plot_base <- function(df.plot,
                      z = "accuracy",
                      contour = TRUE,
                      col.contour = col.clusters,
                      col.pal = rev(pals::brewer.rdbu(10)),
                      z.breaks = acc.breaks.base,
                      z.labels = acc.labels.base,
                      legend.position = "bottom") {
  
  # Onset lines
  lines.v.pos <- 330*seq(round(min(df.plot$t_test)/333), round(max(df.plot$t_test)/333), 1)
  lines.v.pos <- lines.v.pos[lines.v.pos!=0]
  lines.h.pos <- 330*seq(round(min(df.plot$t_train)/333), round(max(df.plot$t_train)/333), 1)
  lines.h.pos <- lines.h.pos[lines.h.pos!=0]

  # Plot
  g <- df.plot %>% 
    ggplot(aes(x = t_test, y = t_train)) +
    # --- heatmap
    geom_tile(aes_string(fill = z)) +
    # --- sounds onsets
    geom_vline(xintercept = lines.v.pos, linetype = 2, color = col.lines) +
    geom_hline(yintercept = lines.h.pos, linetype = 2, color = col.lines) +
    geom_vline(xintercept = 0, color = col.lines) +
    geom_hline(yintercept = 0, color = col.lines) +
    
    scale_x_continuous(breaks = scales::pretty_breaks(12),
                       expand = c(0,0)) +
    scale_y_continuous(breaks = scales::pretty_breaks(4),
                       expand = c(0,0)) +
    scale_fill_stepsn(colors = col.pal,
                      breaks = z.breaks,
                      labels = z.labels,
                      limits = c(min(z.breaks), max(z.breaks)),
                      oob = scales::squish) +
    
    labs(x = "Test time (ms)", y = "Train time (ms)") +
  
    coord_equal(clip = "off") +
    
    theme(# Text size
          plot.title = element_text(size = 12, hjust = 0.5, face = "bold"),
          plot.subtitle = element_text(size = 9, hjust = 0.5),
          axis.title = element_text(size = 9),
          legend.title = element_text(size = 9),
          strip.text = element_text(size = 9),
          axis.text = element_text(size = 5),
          legend.text = element_text(size = 5.5),
          # Time ticks
          axis.ticks = element_line(color = "black"),
          axis.ticks.length = unit(0.1, "cm"),
          # Legend
          legend.position = legend.position,
          legend.box.margin = margin(0,0,0,0),
          legend.box.spacing = unit(0, "pt"),
          # no grid
          panel.grid = element_blank(),
          # don't clip facet titles
          strip.clip = "off")
  
  if (legend.position %in% c("bottom","top")) {
    g <- g + 
      guides(fill = guide_colorbar(barwidth = 0.5*length(z.breaks-1), barheight = 0.5,
                                   ticks = F, title.vjust = 1))
  } else {
    g <- g + 
      guides(fill = guide_colorbar(barheight = 0.5*length(z.breaks-1), barwidth = 0.5,
                                   ticks = F, title.vjust = 1))
  }
  
  # --- significant clusters
  if (contour) {
    g <- g +
      geom_contour(aes(z = signif), size = 0.1, color = col.contour)
  }
  
  return(g)
}

Replication & Reordering (empirical null)

Load data

Sounds

tic()

# Load subject-level data (used for clustering) and derive grand average (used for visualization)
df.stim.mean <- tibble()
df.stim.subj <- tibble()
for (manip in c("","_reord")) { 
  
  for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
    tmp <- load_loop_subj(folder = "sounds",
                          condition = paste0(direction, manip))
    
    df.stim.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
    df.stim.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  }
}

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.

Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(row.name)

  # Now:
  data %>% select(all_of(row.name))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.

df.stim.mean %<>% mutate(stim = "sounds")
df.stim.subj %<>% mutate(stim = "sounds") %>%
  # --- center decoding accuracy on 0
  mutate(array = map(array, .f=\(x){x-0.25}))

toc() # ~1s

1.253 sec elapsed

Omissions

tic()

df.omit.mean <- tibble()
df.omit.subj <- tibble()
for (manip in c("","_reord")) { 
  
  for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
      
    tmp <- load_loop_subj(folder = "omissions",
                          condition = paste0(direction, manip))
    
    df.omit.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
    df.omit.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  }
}

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.

Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(row.name)

  # Now:
  data %>% select(all_of(row.name))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.

./data/results/11920812IONP/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19750430PNRK/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19800616MRGU/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19810726GDZN/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19810918SLBR/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19821223KRHR/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19830114RFTM/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19851130EIFI/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19861231ANSR/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19871229CRBC/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19880328AGSG/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19891222GBHL/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19900606KTAD/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19901026KRKE/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19910823SSLD/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19920804CRLE/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19921111BRHC/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930118IMSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930423EEHB/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930524GNTA/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930621ATLI/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19930630MNSU/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19940403HIEC/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19940601IGSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19940930SZDT/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19950212BTKC/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19950326IIQI/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19950604LLZM/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19950905MCSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19960304SBPE/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19960418GBSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19960708HLHY/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!
./data/results/19961118BRSH/omissions/cv_rd_to_rd_reord_scores.npy  does not exist !!!

df.omit.mean %<>% mutate(stim = "omissions") %>% filter(abs(t_test) <= 670)
df.omit.subj %<>% mutate(stim = "omissions") %>%
  # --- center decoding accuracy on 0
  mutate(array = map(array, .f=\(x){x-0.25})) %>% 
  # --- crop subject-level data for cluster analysis
  mutate(array = map(array, .f=\(x){x[71:104,4:138,]}))

toc() # ~ 2s

3.59 sec elapsed

Correlations with entropy

## Prepare input data
df.tmp <- bind_rows(df.stim.subj, df.omit.subj) %>%
  calculate_diff(f.map = map2, values_to = "array") %>%
  mutate(regularity = str_remove(direction, "rd_to_"),
         regularity = factor(regularity, levels = c("rd","mm","mp","or"))) %>% 
  select(-direction) %>% 
  arrange(as.numeric(regularity))
df.tmp

# A tibble: 24 × 4
   stim      manip     array                 regularity
   <chr>     <chr>     <list>                <fct>     
 1 sounds    original  <dbl [33 × 135 × 33]> rd        
 2 sounds    reordered <dbl [33 × 135 × 33]> rd        
 3 sounds    diff      <dbl [33 × 135 × 33]> rd        
 4 omissions original  <dbl [34 × 135 × 33]> rd        
 5 omissions reordered <dbl [34 × 135 × 33]> rd        
 6 omissions diff      <dbl [34 × 135 × 33]> rd        
 7 sounds    original  <dbl [33 × 135 × 33]> mm        
 8 sounds    reordered <dbl [33 × 135 × 33]> mm        
 9 sounds    diff      <dbl [33 × 135 × 33]> mm        
10 omissions original  <dbl [34 × 135 × 33]> mm        
# ℹ 14 more rows

tic()
df.corr.mean <- tibble()
df.corr.subj <- tibble() 

## Loop through condition "manip" and stimulus type "stim
for (mm in c("original","reordered","diff")) { 

  for (ss in c("sounds", "omissions")) { #
    cat("\n\nProcessing ", str_to_upper(mm), " / ", str_to_upper(ss), "...\n", sep="")
    
    # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
    for (idx.subj in 1:n.subj) {
      
      # --- bind time-generalization matrices for the 3 non-random conditions into a single 3D array
      # (add accuracy data from the random condition if we are processing sounds)
      if (ss == "sounds") {
        idx.lines = 1:4
      } else if (ss == "omissions") {
        idx.lines = 2:4
      }
      arr.tmp <- correlations_extract3D(data = filter(df.tmp, manip == mm, stim == ss),
                                        idx.lines = idx.lines, idx.subj)
      
      # --- calculate correlation with entropy level coded as 0,1,2(,3), obtaining a 2D array
      if (ss == "sounds") {
        arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:3)
      } else if (ss == "omissions") {
        arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 1:3)
      }
      
      # --- append to already processed subjects in a 3D array
      if (idx.subj == 1) {
        arr.corr <- arr.corr.tmp
      } else {
        arr.corr %<>% abind(arr.corr.tmp, along = 3)
      }
    }
    
    # Postprocessing
    arr.corr %<>% 
      correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
      atanh # apply Fisher's transformation
    
    # --- store participant-level results in a dataframe
    df.corr.subj %<>% bind_rows(tibble(manip = mm, stim = ss, array = list(arr.corr)))
    
    # --- calculate groupe average and store in a dataframe
    df.corr.mean %<>% bind_rows(
      arr.corr %>% 
        apply(c(1,2), mean) %>% # average across participants
        tanh %>% # convert back to correlation coefficients  
        matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>% 
        mutate(manip = mm, stim = ss)
    )
  }
}

Processing ORIGINAL / SOUNDS...


Processing ORIGINAL / OMISSIONS...


Processing REORDERED / SOUNDS...


Processing REORDERED / OMISSIONS...


Processing DIFF / SOUNDS...


Processing DIFF / OMISSIONS...

toc() # ~20 s

33.512 sec elapsed

Clusters stats

Accuracy

tic()

df.tmp.subj <- bind_rows(df.stim.subj, df.omit.subj)
df.tmp.mean <- bind_rows(df.stim.mean, df.omit.mean)
  
# Cluster analysis
df.clusters.0 <- df.tmp.subj %>% 
  calculate_diff(f.map = map2, values_to = "array") %>% 
  cluster_loop_cond(n = n_permutations) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))

Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(values_to)

  # Now:
  data %>% select(all_of(values_to))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.

Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(row.name)

  # Now:
  data %>% select(all_of(row.name))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.

# Quick check of clusters with a mask plot
df.clusters.0 %>% 
  mutate(direction = (direction %>% str_replace("rd_to_", "") %>% str_to_upper)) %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid(direction ~ manip) +
  geom_tile() +
  coord_equal()

# Aggregate all data and clean
df.plot.0 <- df.tmp.mean %>%
  calculate_diff(f.map = map2_dbl, values_to = "accuracy") %>% 
  # --- append cluster significance
  left_join(df.clusters.0) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
                                diff = "ORIGINAL – REORDERED"))) %>% 
  mutate(direction = toupper(sub("rd_to_", "", direction)),
         direction = factor(direction, levels = rev(c("RD","MM","MP","OR"))))

Joining with `by = join_by(t_train, t_test, direction, stim, manip)`

toc() # ~3.5 minutes

23.236 sec elapsed

Correlations

tic()

# Cluster analysis
df.clusters.corr <- df.corr.subj %>% 
  cluster_loop_cond(n = n_permutations)

# Quick check of clusters with a mask plot
df.clusters.corr %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid( ~ manip) +
  geom_tile() +
  coord_equal()

# --- append cluster significance
df.plot.0.corr <- df.corr.mean %>% 
  left_join(df.clusters.corr) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
                                diff = "ORIGINAL – REORDERED"))) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))

Joining with `by = join_by(t_train, t_test, manip, stim)`

df.plot.0.corr %>% summary

    t_train          t_test           r                             manip     
 Min.   :  0.0   Min.   :-670   Min.   :-0.91933   ORIGINAL            :9045  
 1st Qu.: 80.0   1st Qu.:-340   1st Qu.:-0.21316   REORDERED           :9045  
 Median :160.0   Median :   0   Median :-0.01320   ORIGINAL – REORDERED:9045  
 Mean   :162.5   Mean   :   0   Mean   :-0.01764                              
 3rd Qu.:250.0   3rd Qu.: 340   3rd Qu.: 0.18399                              
 Max.   :330.0   Max.   : 670   Max.   : 0.86931                              
     stim               signif      
 Length:27135       Min.   :0.0000  
 Class :character   1st Qu.:0.0000  
 Mode  :character   Median :0.0000  
                    Mean   :0.1094  
                    3rd Qu.:0.0000  
                    Max.   :1.0000  

toc() # ~ 1 min

3.179 sec elapsed

Theoretical null

Load confusion matrices & calculate

tic()

null_accuracy <- function(mat.confusion, mat.transition) {
  # Calculate accuracy under the null hypothesis from the transition and confusion matrices mat.T and mat.C
  # The 3/4 factor is used to compensate the unbiased estimator used in cov() by default (see ?cov)
  return(0.25 + sum(diag(3/4*cov(t(mat.transition), t(mat.confusion)))))
}

df.maths.subj <- tibble()
df.maths.mean <- tibble()
for (entropy in c("RD","MM","MP","OR")) {

  array.scores <- array(0.25, dim = c(33,136,n.subj))
  
  i <- 1
  for (s in list.subj) {
    
    # --- load the confusion matrix
    tmp.file <- os$path$join(path.root, s, "rd_to_rd_confmats.npz") %>%
      np$load()
    
    # --- crop around S0 (decoded stimulus)
    mat.C <- tmp.file$f[["arr_0"]][,71:103,71:104,,] %>%
      apply(c(2,3,4,5), mean)
    
    # NOTE: because the transition matrix is the same for all subjects, theoretical accuracy can be calculated from the group average of the confusion matrix (faster)
    # --- calculate theoretical accuracy
    #   --- for S-2
    array.scores[1:33,1:34,i] <- mat.C %>%
      apply(c(1,2), null_accuracy, mat.transition = mat.T[[entropy]]%^%2)
    #   --- for S-1
    array.scores[1:33,35:68,i] <- mat.C %>%
      apply(c(1,2), null_accuracy, mat.transition = mat.T[[entropy]]%^%1)
    #   --- for S
    array.scores[1:33,69:102,i] <- mat.C %>%
      apply(c(1,2), \(x)mean(diag(x)))
    #   --- for S+1
    array.scores[1:33,103:136,i] <- mat.C %>%
      apply(c(1,2), null_accuracy, mat.transition = t(mat.T[[entropy]]))
    i <- i+1
  }
  
  array.scores <- array.scores[,2:136,]
  
  # Append to other entropy levels in dataframe
  # --- subject-level data
  df.maths.subj %<>% bind_rows(tibble(array = list(array.scores), entropy = entropy))
  # --- group-averaged data
  df.maths.mean %<>% bind_rows(array.scores %>% 
      apply(c(1,2), mean) %>%
      matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "accuracy") %>% 
    mutate(entropy = entropy)
  )
}

toc() # ~ 1min

67.387 sec elapsed

Append original data & calculate diff

df.maths.all <- 
# --- aggregate original, theoretical and difference between the 2
  bind_rows(df.stim.subj %>%
              # --- select original analysis
              filter(manip == "") %>% mutate(manip = "original") %>% 
              # --- crop to the same time-window than for theoretical analysis
              mutate(entropy = toupper(str_replace(direction,"rd_to_",""))),
            df.maths.subj %>%
              mutate(manip = "maths") %>% 
              mutate(array = map(array, .f=\(x){x-0.25}))) %>% 
  select(manip, entropy, array) %>% 
  pivot_wider(names_from = manip, values_from = "array") %>% 
  mutate(diff = map2(original, maths, .f=\(x,y){x-y})) %>% 
  pivot_longer(cols = c(original, maths, diff), names_to = "manip", values_to = "array")

Correlations with entropy

tic()

df.corr.maths.mean <- tibble()
df.corr.maths.subj <- tibble() 

for (mm in c("original","maths","diff")) { 
  # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
  for (idx.subj in 1:n.subj) {
    
    # --- bind time-generalization matrices of the 4 conditions into a single 3D array
    arr.tmp <- correlations_extract3D(data = df.maths.all %>% filter(manip == mm),
                                      idx.lines = 1:4, idx.subj)
    
    # --- calculate correlation with entropy level coded as 0,1,2,3, obtaining a 2D array
    arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:3)
    
    # --- append to already processed subjects in a 3D array
    if (idx.subj == 1) {
      arr.corr <- arr.corr.tmp
    } else {
      arr.corr %<>% abind(arr.corr.tmp, along = 3)
    }
  }
  
  # Postprocessing
  arr.corr %<>% 
    correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
    atanh # apply Fisher's transformation
    
  # --- store participant-level results in a dataframe of arrays
  df.corr.maths.subj %<>% bind_rows(tibble(manip = mm, array = list(arr.corr)))
  
  # --- calculate group average and store in a dataframe
  df.corr.maths.mean %<>% bind_rows(
    arr.corr %>% 
      apply(c(1,2), mean) %>% # average across participants
      tanh %>% # convert back to correlation coefficients  
      matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>% 
      mutate(manip = mm)
  )
}

toc() # ~ 20s

16.723 sec elapsed

Clusters stats

Accuracy

tic()

df.clusters.maths <- df.maths.all %>% 
  cluster_loop_cond(n = n_permutations) %>%
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))

toc() # ~ 2 min

2.065 sec elapsed

# Quick check of clusters with a mask plot
df.clusters.maths %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid(entropy ~ manip) +
  geom_tile() +
  coord_equal()

Correlations

tic()
df.clusters.corr.maths <- df.corr.maths.subj %>% 
  cluster_loop_cond(n = n_permutations)
toc() # ~ 20s

0.566 sec elapsed

# Quick check of clusters with a mask plot
df.clusters.corr.maths %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid( ~ manip) +
  geom_tile() +
  coord_equal()

# --- append cluster significance
df.plot.maths.corr <- df.corr.maths.mean %>% 
  left_join(df.clusters.corr.maths) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", maths = "THEORETICAL NULL",
                                diff = "ORIGINAL – THEORETICAL NULL"))) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))

Joining with `by = join_by(t_train, t_test, manip)`

df.plot.maths.corr %>% summary

    t_train        t_test           r           
 Min.   :  0   Min.   :-670   Min.   :-0.98948  
 1st Qu.: 80   1st Qu.:-340   1st Qu.:-0.23120  
 Median :160   Median :   0   Median : 0.00000  
 Mean   :160   Mean   :   0   Mean   :-0.03176  
 3rd Qu.:240   3rd Qu.: 340   3rd Qu.: 0.16949  
 Max.   :320   Max.   : 670   Max.   : 0.95706  
                         manip          signif      
 ORIGINAL                   :4455   Min.   :0.0000  
 THEORETICAL NULL           :4455   1st Qu.:0.0000  
 ORIGINAL – THEORETICAL NULL:4455   Median :0.0000  
                                    Mean   :0.1631  
                                    3rd Qu.:0.0000  
                                    Max.   :1.0000  

Prediction of the most likely (“simple prediction”)

Load data

tic()

# Load subject-level data and derive grand average
df.stim.ml.mean <- tibble()
df.stim.ml.subj <- tibble()
for (manip in c("","_reord")) { 
  
  for (direction in c("rd_to_mm", "rd_to_mp", "rd_to_or")) {
    
    tmp <- load_loop_subj(folder = "sounds",
                          condition = paste0(direction, manip, "_sp"))
    
    df.stim.ml.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
    df.stim.ml.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  }
}

df.stim.ml.mean %<>% mutate(stim = "sounds")
df.stim.ml.subj %<>% mutate(stim = "sounds") %>%
  # --- center decoding accuracy on 0
  mutate(array = map(array, .f=\(x){x-0.25}))

toc() # ~ 1s

0.891 sec elapsed

df.stim.ml.mean %>% mutate_if(is.character, as.factor) %>% summary

    t_train        t_test        accuracy         manip          direction   
 Min.   :  0   Min.   :-670   Min.   :0.2040         :13365   rd_to_mm:8910  
 1st Qu.: 80   1st Qu.:-340   1st Qu.:0.2465   _reord:13365   rd_to_mp:8910  
 Median :160   Median :   0   Median :0.2498                  rd_to_or:8910  
 Mean   :160   Mean   :   0   Mean   :0.2500                                 
 3rd Qu.:240   3rd Qu.: 340   3rd Qu.:0.2532                                 
 Max.   :320   Max.   : 670   Max.   :0.3095                                 
     stim      
 sounds:26730  
               
               
               
               
               

Correlations with entropy

tic()

## Prepare input data
df.tmp <- df.stim.ml.subj %>% 
  calculate_diff(f.map = map2, values_to = "array") %>%
  mutate(entropy = str_remove(direction, "rd_to_"),
         entropy = factor(entropy, levels = c("rd","mm","mp","or"))) %>% 
  select(-direction) %>% 
  arrange(as.numeric(entropy))


## Calculate correlations
df.corr.ml.mean <- tibble()
df.corr.ml.subj <- tibble() 

# --- Loop through condition "manip" and stimulus type "stim"
for (mm in c("original","reordered","diff")) { 

  for (ss in c("sounds")) {
    cat("\n\nProcessing ", str_to_upper(mm), " / ", str_to_upper(ss), "...\n", sep="")

    
    # --- Loop through subjects to derive a 3D array of participant-level time-generalized correlations with entropy, from a 3D array of participant-level accuracy matrix
    for (idx.subj in 1:n.subj) {
      
      # --- bind time-generalization matrices for the 3 non-random conditions into a single 3D array
      arr.tmp <- correlations_extract3D(data = filter(df.tmp, manip == mm, stim == ss),
                                        idx.lines = 1:3, idx.subj)
      
      # --- calculate correlation with entropy level coded as 0,1,2, obtaining a 2D array
      arr.corr.tmp <- apply(arr.tmp, 1:2, correlations_vecs, vec2 = 0:2)
      
      # --- append to already processed subjects in a 3D array
      if (idx.subj == 1) {
        arr.corr <- arr.corr.tmp
      } else {
        arr.corr %<>% abind(arr.corr.tmp, along = 3)
      }
    }
    
    # Post-processing
    arr.corr <- arr.corr %>% 
      correlations_fixval %>% # replace NAs by 0 and extreme values by +/- .99
      atanh # apply Fisher's transformation
    
    # --- store participant-level results in a dataframe
    df.corr.ml.subj %<>% bind_rows(tibble(manip = mm, stim = ss, array = list(arr.corr)))
    
    # --- calculate groupe average and store in a dataframe
    df.corr.ml.mean %<>% bind_rows(
      arr.corr %>% 
        apply(c(1,2), mean) %>% # average across participants
        tanh %>% # convert back to correlation coefficients  
        matrix_to_df2(row.name = "t_train", col.name = "t_test", val.name = "r") %>% 
        mutate(manip = mm, stim = ss)
    )
  }
}

Processing ORIGINAL / SOUNDS...


Processing REORDERED / SOUNDS...


Processing DIFF / SOUNDS...

toc() # ~ 15s

11.529 sec elapsed

Cluster stats

Accuracy

tic()

## Cluster analysis with subject-level data
# --- calculate original-reordered difference
tmp.diff <- df.stim.ml.subj %>% 
  calculate_diff(f.map = map2, values_to = "array")

# --- combine & run cluster analysis
df.clusters.ml <- tmp.diff %>% #bind_rows(tmp.diff, tmp.ave) %>% 
  cluster_loop_cond(n = n_permutations) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))

# Quick check of clusters with a mask plot
df.clusters.ml %>% 
  mutate(direction = (direction %>% str_replace("rd_to_", "") %>% str_to_upper)) %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid(direction ~ manip) +
  geom_tile() +
  coord_equal()

## Grand average data
# --- calculate original-reordered difference
tmp.diff <- df.stim.ml.mean %>%
  calculate_diff(f.map = map2_dbl, values_to = "accuracy")
# --- calculate average of difference across entropy levels
tmp.ave <- tmp.diff %>% 
  filter(manip == "diff") %>% 
  group_by(t_train, t_test, manip, stim) %>% 
  summarise(accuracy = mean(accuracy, na.rm=T)) %>% 
  mutate(direction = "average of\n(MM,MP,OR)")

`summarise()` has grouped output by 't_train', 't_test', 'manip'. You can
override using the `.groups` argument.

# --- combine
df.plot.ml <- tmp.diff #bind_rows(tmp.diff, tmp.ave)

##  Aggregate all data and clean
df.plot.ml %<>%
  # --- append cluster significance
  left_join(df.clusters.ml) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
                                diff = "ORIGINAL – REORDERED"))) %>%
  mutate(direction = toupper(sub("rd_to_", "", direction))) %>% 
  mutate(direction = factor(direction, levels = rev(c("AVERAGE OF\n(MM,MP,OR)","MM","MP","OR"))))

Joining with `by = join_by(t_train, t_test, direction, stim, manip)`

toc() # ~1.5min

2.544 sec elapsed

Correlations

tic()

df.clusters.corr.ml <- df.corr.ml.subj %>% 
  cluster_loop_cond(n = n_permutations)

toc() # ~ 30s

0.583 sec elapsed

# Quick check of clusters with a mask plot
df.clusters.corr.ml %>% 
  ggplot(aes(x = t_test, y = t_train, fill = signif)) +
  facet_grid( ~ manip) +
  geom_tile() +
  coord_equal()

# --- append cluster significance
df.plot.ml.corr <- df.corr.ml.mean %>% 
  filter(between(t_test,35,68)) %>% 
  left_join(df.clusters.corr.ml) %>%
  mutate(manip = (recode_factor(manip, original = "ORIGINAL", reordered = "REORDERED",
                                diff = "ORIGINAL – REORDERED"))) %>% 
  # --- convert indices to time
  mutate(t_test = (dt*(t_test-1)-time_offset)) %>% 
  mutate(t_train = (dt*(t_train-1)))

Joining with `by = join_by(t_train, t_test, manip, stim)`

df.plot.ml.corr %>% summary

    t_train        t_test           r                             manip     
 Min.   :  0   Min.   :-330   Min.   :-0.75166   ORIGINAL            :1122  
 1st Qu.: 80   1st Qu.:-250   1st Qu.:-0.24421   REORDERED           :1122  
 Median :160   Median :-165   Median :-0.05519   ORIGINAL – REORDERED:1122  
 Mean   :160   Mean   :-165   Mean   :-0.05098                              
 3rd Qu.:240   3rd Qu.: -80   3rd Qu.: 0.13428                              
 Max.   :320   Max.   :   0   Max.   : 0.77497                              
     stim               signif         
 Length:3366        Min.   :0.0000000  
 Class :character   1st Qu.:0.0000000  
 Mode  :character   Median :0.0000000  
                    Mean   :0.0008913  
                    3rd Qu.:0.0000000  
                    Max.   :1.0000000  

Figures

Methodological plots

Transition matrices (fig. 1a)

# Convert & combine transition matrices into a dataframe
df.plot <- tibble()
for (condition in c("RD","MM","MP","OR")) {
  df.plot %<>% bind_rows(
    mat.T[[condition]] %>% 
      matrix_to_df2(row.name = "from", col.name = "to", val.name = "p") %>% 
      mutate(condition = condition)
  )  
}
df.plot %<>%
  mutate(condition = factor(condition, names(mat.T)))

df.plot.pitch <- bind_rows(tibble(x = 0.1, y = 1:4, pitch = 1:4),
                           tibble(y = 0.1, x = 1:4, pitch = 1:4))
# Plot
df.plot %>% 
  ggplot(aes(x = to, y = from)) +
  facet_wrap(~ condition, nrow = 1, strip.position = "bottom") +
  # --- heatmap
  geom_tile(aes(fill = p), color = "black") +
  # --- display percentages
  geom_text(data = . %>% filter(p<=.25),
            aes(label = paste0(100*p,"%")),
            color = "black", size = 2) +
  geom_text(data = . %>% filter(p>.25),
            aes(label = paste0(100*p,"%")),
            color = "white", size = 2) +
  
  # --- display color-coded pitch in axes
  geom_point(data = df.plot.pitch,
             aes(x, y, color = pitch), size = 4, show.legend = FALSE) +
  geom_text(data = df.plot.pitch,
             aes(x, y, label = LETTERS[pitch]), color = "white", size = 2, show.legend = FALSE) +
  
  scale_color_viridis_c(end = 0.8) +
  scale_x_continuous(name = "to\n", position = "top") +
  scale_y_reverse(name = "from\n") +
  scale_fill_gradient(name = "transition probability",
                      low = "white", high = "black",
                      labels = ~str_replace(., "^0.","."),
                      limits = c(0,1)) +
  guides(fill = guide_colorbar(barheight = 0.5,
                               ticks = F, title.vjust = 1)) +
  coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
              expand = F, clip = "off") +
  theme(legend.position = "bottom",
        legend.box.margin = margin(t = -20),
        # text size
        axis.title = element_text(size = 9),
        legend.title = element_text(size = 9),
        strip.text = element_text(size = 12),
        # remove axis text
        axis.text = element_blank(),
        # Positions & margins
        axis.title.x = element_text(hjust = 0.08),
        panel.spacing.x = unit(0.8,"cm")) -> g.transitions

g.transitions

Reordering (fig. 1b)

# Generate Markov sequences
n.samples <- 12
stims <- list(RD = c(), OR = c())
set.seed(128764) #20476, 99713
for (entropy in c("RD","OR")) {
  stims[[entropy]] <- vector("numeric", n.samples)
  stims[[entropy]][1] <- sample(1:4, 1)
  for (i in 2:n.samples) {
    stims[[entropy]][i] <- sample(1:4, 1, prob = mat.T[[entropy]][stims[[entropy]][i-1],])
  }
}

# Add global sequence indices, class-specific indices and offset according to pitch
df.plot <- stims %>% as.tibble() %>% 
  mutate(idx = as.double(1:n.samples)) %>% 
  pivot_longer(-idx, names_to = "entropy", values_to = "class") %>% 
  mutate(entropy = factor(entropy, levels = c("OR","RD"))) %>%
  group_by(entropy,class) %>% mutate(rank = rank(idx)) %>% 
  mutate(y = as.double(entropy) + class/16-10/64)

pitch.y <- sort(unique(df.plot$y))

# Plot
df.plot %>% 
  ggplot(aes(x = idx, y = y, color = class)) +
  
  # --- pitch lines
  geom_hline(yintercept = pitch.y, color = "grey90") +
  # --- pitch labels
  annotate(geom = "text", label = "pitch",
           x = 12.5, y = c(mean(head(pitch.y,4)), mean(tail(pitch.y,4))),
           hjust = 0.5, vjust = 2.5, size = 2.5, color = "grey50", angle = 90) +
  
  annotate(geom = "text", label = "low",
           x = 12.5, y = c(min(head(pitch.y,4)), min(tail(pitch.y,4))),
           hjust = 0.74, vjust = 1.2, size = 2, color = "grey50", angle = 90) +
  
  annotate(geom = "text", label = "high",
           x = 12.5, y = c(max(head(pitch.y,4)), max(tail(pitch.y,4))),
           hjust = 0.25, vjust = 1.2, size = 2, color = "grey50", angle = 90) +
  
  # --- reordering trajectories
  ggh4x::geom_pointpath(data = df.plot %>% arrange(desc(entropy)),
               aes(x = idx, y = y, group = interaction(class,rank)),
               mult = 0.4, color = "black", linewidth = 0.3,
               arrow = arrow(type = "closed", angle = 25, length=unit(.15, 'cm'))) +
  # --- stimuli
  geom_point(size = 4, show.legend = F) +
  # --- stimuli ranks
  geom_text(aes(label = rank), size = 2, color = "white") +
  
  scale_x_continuous(name = "trial number", breaks = 1:n.samples) +
  scale_y_continuous(breaks = c(1,2), labels = c("OR","RD"), position = "left") +
  scale_color_viridis_c(end = 0.8) +
  coord_cartesian(clip = "off") +
  
  theme(axis.title.y = element_blank(),
        axis.text.y = element_text(size = 10, color = "black", margin = margin(r=-2)),
        axis.title.x = element_text(size = 9, color = "black"),
        panel.grid = element_blank(),
        panel.grid.major.x = element_line(color = "grey80", linewidth = 0.2)) -> g.reordering

g.reordering

Prediction types (fig. 2a)

Presented vs. most likely

# Generate Markov sequences
n.samples <- 6
stims <- mat.T
set.seed(12871) 
for (entropy in names(stims)) {
  stims[[entropy]] <- vector("numeric", n.samples)
  stims[[entropy]][1] <- 4#sample(1:4, 1)
  for (i in 2:n.samples) {
    stims[[entropy]][i] <- sample(1:4, 1, prob = mat.T[[entropy]][stims[[entropy]][i-1],])
  }
}

# Add sequence index & predictions
df.plot <- stims %>% as.tibble() %>% 
  mutate(idx = as.double(1:n.samples)) %>% 
  pivot_longer(-idx, names_to = "entropy", values_to = "class") %>% 
  mutate(entropy = factor(entropy, levels = rev(names(stims)))) %>% 
  group_by(entropy) %>% mutate(actual = lead(class),
                               mostlikely = ifelse(class==1, 4, class-1)) %>% 
  pivot_longer(c(actual,mostlikely), names_to = "status", values_to = "class.pred") %>% 
  filter(entropy != "RD")

# Plot
df.plot %>% 
  ggplot(aes(x = idx, y = class)) +
  facet_wrap(~ entropy, ncol = 1, strip.position = "right") +
  
  geom_segment(data = df.plot %>% filter(status == "mostlikely", idx != n.samples),
               aes(x = idx, xend = idx+1, y = class, yend = class.pred, linetype = as.factor(sign(class.pred-class))),
               #arrow = arrow(type = "closed", length = unit(6,"pt")),
               color = "grey50", show.legend = F) +
  
  geom_point(aes(fill = class), size = 3, shape = 21, color = "white", show.legend = T) +
  
  geom_point(data = df.plot %>% filter(status == "mostlikely", idx != n.samples),
             aes(x = idx+1, y = class.pred, color = class.pred),
             size = 2.5, shape = 8) +
  
  scale_x_continuous(name = "trial number", breaks = 1:n.samples) +
  scale_y_continuous(name = "") +
  scale_color_viridis_c(end = 0.8, guide = "legend", name = "", labels = c("","","","most likely")) +
  scale_fill_viridis_c(end = 0.8, guide = "legend", name = "", labels = c("","","","actually presented")) +
  
  guides(linetype = F) +
  coord_cartesian(clip = "off") +
  expand_limits(x = c(0.9,6.1), y = c(0.3,4.7)) +

  guides(#color = guide_legend(order = 1, override.aes = list(color = "grey50")),
         color = guide_legend(order = 1),
         fill = guide_legend(order = 2)) +
  theme(plot.caption = element_text(size = 6),
        axis.text.y = element_blank(),
        axis.title.x = element_text(size = 9, color = "black"),
        strip.text = element_text(size = 12),
        panel.spacing = unit(0.2,"cm"),
        panel.grid = element_blank(),
        panel.border = element_rect(fill=NA),
        legend.position = "top",
        legend.box = "vertical",
        legend.box.just = "left",
        legend.margin = margin(b = -10),
        legend.text = element_text(margin = margin(r = -12, l=0,  unit = "pt"))
        ) -> g.mostlikely

g.mostlikely

Figure 1 & Supplementary Figure 2

Initialize

stims <- list(RD = c(), OR = c())
g.base <- list(sound = c(), omission = c())
g.diff <- list(sound = c(), omission = c())

# Create ad hoc dataframes for onset' labels
df.labels.onset <- list(
  sounds = tibble(label = "0 = sound onset",
                manip = factor(c("REPLICATION\n", "EMPIRICAL NULL\n",
                                 "REPLICATION – EMPIRICAL NULL\n")),
                direction = factor("OR", levels = names(stims))),
  omissions = tibble(label = "0 = omission onset",
                manip = factor(c("REPLICATION\n", "EMPIRICAL NULL\n",
                                 "REPLICATION – EMPIRICAL NULL\n")),
                direction = factor("OR", levels = names(stims)))
)
  
# Function to add onset' labels to a plot
add_onset <- function(g, df.labels) {
  g <- g + geom_text(data = df.labels %>% filter(manip %in% unique(g$data$manip)),
                     aes(label = label, x = 0, y = 350),
                     hjust = 0, vjust = -0.5, size = 3, fontface = "italic")
  return(g)
}

Original & empirical null: accuracy plots

for (s in c("sounds","omissions")) { #
  # ORIGINAL & REORDERED
  df.plot.0 %>% 
    filter(stim == s, manip != "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip, 
                                 ORIGINAL = "REPLICATION\n", 
                                 REORDERED = "EMPIRICAL NULL\n")) %>%
    plot_base(col.pal = col.pal.base,
              z.breaks = acc.breaks.base,
              z.labels = acc.labels.base) +
    facet_grid(direction ~ manip) -> g
  
  g.base[["acc"]][[s]] <- g %>% add_onset(df.labels.onset[[s]])
  
  # DIFFERENCE
  g.diff[["acc"]][[s]] <- df.plot.0 %>%
    filter(stim == s, manip == "ORIGINAL – REORDERED") %>%
    mutate(manip = recode_factor(manip,
                                 `ORIGINAL – REORDERED` = "REPLICATION – EMPIRICAL NULL\n")) %>%
    plot_base(col.pal = col.pal.diff,
              z.breaks = acc.breaks.diff,
              z.labels = acc.labels.diff) +
    facet_grid(direction ~ manip) +
    theme(axis.title.y = element_blank()) -> g
  
  g.diff[["acc"]][[s]] <- g %>% add_onset(df.labels.onset[[s]])
}

# Assemble Figure 1c
g.fig1.1.acc <- (
  g.base[["acc"]][["sounds"]] + 
    (g.diff[["acc"]][["sounds"]] + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1)) &
    theme(strip.text.x = element_text(face = "bold"))
)

g.fig1.1.acc

# Assemble Supplementary Figure 1a
g.supfig2.acc <- (
  g.base[["acc"]][["omissions"]] + 
    (g.diff[["acc"]][["omissions"]] + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1)) &
    theme(strip.text.x = element_text(face = "bold"))
)

g.supfig2.acc

Original & empirical null: correlations plots

for (s in c("sounds", "omissions")) {
  # ORIGINAL & REORDERED
  df.plot.0.corr %>% 
    filter(stim == s, manip != "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip, 
                                 ORIGINAL = "REPLICATION\n", #  (original dataset)\n
                                 REORDERED = "EMPIRICAL NULL\n")) %>%  # (reordered RD)\n
    mutate(direction = "") %>% 
    plot_base(col.pal = col.pal.corr,
              z = "r",
              z.breaks = corr.breaks,
              z.labels = corr.labels) +
    labs(fill = "Correlation") +
    facet_grid(direction ~ manip) -> g
  
  g.base[["corr"]][[s]] <- g #%>% add_onset(df.labels.onset[[s]] %>% mutate(direction = ""))
  
  # DIFFERENCE
  df.plot.0.corr %>%
    filter(stim == s, manip == "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip,
                                 `ORIGINAL – REORDERED` = "REPLICATION – EMPIRICAL NULL\n")) %>% 
    mutate(direction = "") %>% 
    plot_base(col.pal = col.pal.corr,
              z = "r",
              z.breaks = corr.breaks,
              z.labels = corr.labels) +
    facet_grid(. ~ manip) +
    labs(fill = "Correlation") +
    theme(axis.title.y = element_blank()) -> g
  
  g.diff[["corr"]][[s]] <- g #%>% add_onset(df.labels.onset[[s]] %>% mutate(direction = ""))
}

# Assemble Figure 1d
g.fig1.1.corr <- (
  g.base[["corr"]][["sounds"]] + 
    (g.diff[["corr"]][["sounds"]] + plot_layout(tag_level = "new")) + 
    plot_layout(widths = c(2, 1), guides = "collect") &
    theme(strip.text = element_blank(),
          legend.position = "bottom",
          legend.margin = margin(t = -20, l = 60))
)

g.fig1.1.corr

# Assemble Supplementary Figure 1a
g.supfig2.corr <- (
  g.base[["corr"]][["omissions"]] + 
    (g.diff[["corr"]][["omissions"]] + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1)) &
    theme(strip.text.x = element_text(face = "bold"))
)

g.supfig2.corr

Theoretical null: accuracy plots

## Aggregate all data and clean

df.plot.maths <- bind_rows( 
  # --- make datasets compatible
  df.maths.mean %>% mutate(manip = "maths") %>% 
    mutate(t_train = dt*(t_train-1) - 0,
           t_test = dt*(t_test-1) - time_offset),
  df.stim.mean %>% 
    filter(manip == "") %>% 
    mutate(manip = "original") %>% 
    mutate(entropy = toupper(sub("rd_to_", "", direction))) %>% 
    # mutate(across(starts_with("t_"), ~(dt*(.-1)-700))) %>% 
    select(-stim, -direction)
) %>% 
  # --- calculate difference
  pivot_wider(names_from = manip, values_from = "accuracy") %>%
  mutate(diff = original - maths) %>%
  pivot_longer(cols = c(original, maths, diff),
               names_to = "manip", values_to = "accuracy") %>% 
  # --- crop time window
  filter(t_train>=0, (t_train-333) < dt) %>% 
  # --- aggregate results of cluster analysis
  left_join(df.clusters.maths) %>% 
  # --- clean labels
  mutate(manip = (recode_factor(manip, original = "REPLICATION\n", maths = "THEORETICAL NULL\n",
                                diff = "REPLICATION – THEORETICAL NULL\n")),
         direction = factor(entropy, levels = rev(c("RD","MM","MP","OR"))))

Joining with `by = join_by(t_train, t_test, entropy, manip)`

## PLOT
# Initialize subplots
g.maths.base <- list()
g.maths.diff <- list()

# --- ORIGINAL & THEORETICAL NULL
df.plot.maths %>% 
  filter(manip != "REPLICATION – THEORETICAL NULL\n") %>% 
  plot_base(col.pal = col.pal.base,
            z.breaks = acc.breaks.base,
            z.labels = acc.labels.base) +
  facet_grid(direction ~ manip) -> g

g.maths.base[["acc"]] <- g %>% add_onset(df.labels.onset[["sounds"]] %>% 
                                           mutate(manip = str_replace(manip, "EMPIRICAL", "THEORETICAL")))

# --- DIFFERENCE
df.plot.maths %>%
  filter(manip == "REPLICATION – THEORETICAL NULL\n") %>% 
  plot_base(col.pal = col.pal.diff,
            z.breaks = acc.breaks.diff,
            z.labels = acc.labels.diff) +
  facet_grid(direction ~ manip) +
  theme(axis.title.y = element_blank()) -> g

g.maths.diff[["acc"]] <- g %>% add_onset(df.labels.onset[["sounds"]] %>% 
                                           mutate(manip = str_replace(manip, "EMPIRICAL", "THEORETICAL")))

g.fig1.2.acc <- (
  g.maths.base[["acc"]] + 
    (g.maths.diff[["acc"]]  + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1)) &
    theme(strip.text.x = element_text(face = "bold"))
  )

g.fig1.2.acc

Theoretical null: correlations plots

df.plot <- df.plot.maths.corr %>% 
  mutate(manip = str_replace(manip, "ORIGINAL", "REPLICATION"),
         manip = paste0(manip, "\n"))

# ORIGINAL & THEORETICAL
df.plot %>% 
  filter(manip != "REPLICATION – THEORETICAL NULL\n") %>% 
  mutate(direction = "") %>% 
  plot_base(col.pal = col.pal.corr,
            z = "r",
            z.breaks = corr.breaks,
            z.labels = corr.labels) +
  labs(fill = "Correlation") +
  facet_grid(direction ~ manip) -> g

g.maths.base[["corr"]] <- g

# DIFFERENCE
df.plot %>%
  filter(manip == "REPLICATION – THEORETICAL NULL\n") %>% 
  mutate(direction = "") %>% 
  plot_base(col.pal = col.pal.corr,
            z = "r",
            z.breaks = corr.breaks,
            z.labels = corr.labels) +
  facet_grid(. ~ manip) +
  labs(fill = "Correlation") +
  theme(axis.title.y = element_blank()) -> g


g.maths.diff[["corr"]] <- g

g.fig1.2.corr <- (
  g.maths.base[["corr"]] + 
    (g.maths.diff[["corr"]] + plot_layout(tag_level = "new")) +
    plot_layout(widths = c(2, 1), guides = "collect") &
    theme(strip.text.x = element_blank(),
          legend.position = "bottom",
          legend.margin = margin(t = -20, l = 60))
  )

g.fig1.2.corr

Assemble Figure 1

((g.transitions + theme(panel.spacing.x = unit(5,"mm"), legend.box.margin = margin(t = -30))) +
   g.reordering + theme(axis.text.x = element_text(margin = margin(t = -15))) +
   plot_layout(widths = c(1.7, 1))) / 
  ((g.fig1.1.acc / g.fig1.1.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) + 
     plot_layout(heights = c(4.7, 1))) / 
  ((g.fig1.2.acc / g.fig1.2.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) + 
     plot_layout(heights = c(4.7, 1))) +
  plot_layout(heights = c(1., 4, 4)) +
  plot_annotation(tag_levels = "a") -> g.fig1

g.fig1

if (save) {
  g.fig1 %>% ggsave(filename = "./figures/fig1_MattersArising.pdf", width = 7, height = 11, device = cairo_pdf)
}

Assemble Supplementary Figure 2 (omissions)

((g.supfig2.acc / g.supfig2.corr & theme(axis.text.x = element_text(angle = 45, hjust = 1))) + 
   plot_layout(heights = c(4, 1))) +
  plot_annotation(tag_levels = "a") -> g.supfig2

g.supfig2

if (save) {
  g.supfig2 %>% ggsave(filename = "./figures/supfig2_MattersArising.pdf", width = 7, height = 5.5, device = cairo_pdf)
  g.supfig2 %>% ggsave(filename = "./figures/supfig2_MattersArising.png", width = 7, height = 5.5, device = ragg::agg_png(), dpi = 600)
}

Figure 2

Initialize

g.base.ml <- list()
g.diff.ml <- list()

Build accuracy plots

for (s in c("sounds")) {
  # ORIGINAL & REORDERED
  df.plot.ml %>% 
    filter(t_test >=-333, t_test <= 0) %>% 
    filter(stim == s, manip != "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip, 
                                 ORIGINAL = "ORIGINAL", 
                                 REORDERED = "EMPIRICAL NULL")) %>%
    plot_base(col.pal = col.pal.base,
              z.breaks = acc.breaks.base,
              z.labels = acc.labels.base,
              legend.position = "right") +
    facet_grid(direction ~ manip) +
    labs(title = paste0("Decoding MOST LIKELY sounds")) +
    theme(plot.title = element_text(hjust = -0.2)) -> g.base.ml[["acc"]][[s]]
  
  # DIFFERENCE
  df.plot.ml %>%
    filter(t_test >=-333, t_test <= 0) %>% 
    filter(!grepl("AVERAGE",direction)) %>% 
    filter(stim == s, manip == "ORIGINAL – REORDERED") %>% 
    mutate(manip = "DIFF") %>% 
    plot_base(col.pal = col.pal.diff,
              z.breaks = acc.breaks.diff,
              z.labels = acc.labels.diff,
              legend.position = "right") +
    facet_grid(direction ~ manip) +
    theme(axis.title.y = element_blank()) -> g.diff.ml[["acc"]][[s]]
}

(g.base.ml[["acc"]][["sounds"]] + (g.diff.ml[["acc"]][["sounds"]] + plot_layout(tag_level = "new")) + plot_layout(widths = c(2, 1))) + 
  theme(plot.tag.position = c(0, 1)) -> g.fig2.acc

g.fig2.acc

Warning: `stat_contour()`: Zero contours were generated

Warning in min(x): no non-missing arguments to min; returning Inf

Warning in max(x): no non-missing arguments to max; returning -Inf

Build correlations plots

for (s in c("sounds")) {
  # ORIGINAL & REORDERED
  df.plot.ml.corr %>% 
    filter(t_test >=-333, t_test <= 0) %>% 
    filter(stim == s, manip != "ORIGINAL – REORDERED") %>% 
    mutate(manip = recode_factor(manip, 
                                 ORIGINAL = "ORIGINAL", 
                                 REORDERED = "EMPIRICAL NULL")) %>%
    mutate(direction = "") %>% 
    plot_base(col.pal = col.pal.corr,
              z = "r",
              z.breaks = corr.breaks,
              z.labels = corr.labels,
              legend.position = "right") +
    labs(fill = "Correlation") +
    facet_grid(direction ~ manip) -> g.base.ml[["corr"]][[s]]
  
  # DIFFERENCE
  df.plot.ml.corr %>%
    filter(t_test >=-333, t_test <= 0) %>% 
    filter(stim == s, manip == "ORIGINAL – REORDERED") %>%
    mutate(manip = "DIFF") %>% 
    mutate(direction = "") %>%
    plot_base(col.pal = col.pal.corr,
              z = "r",
              z.breaks = corr.breaks,
              z.labels = corr.labels,
              legend.position = "right") +
    facet_grid(. ~ manip) +
    labs(fill = "Correlation") +
    theme(axis.title.y = element_blank()) -> g.diff.ml[["corr"]][[s]]
}

(g.base.ml[["corr"]][["sounds"]] + g.diff.ml[["corr"]][["sounds"]] + plot_layout(widths = c(2, 1))) +
  # plot_annotation(tag_levels = list("b.")) &
  theme(plot.tag.position = c(0, 1)) -> g.fig2.corr

g.fig2.corr

Warning: `stat_contour()`: Zero contours were generated

Warning in min(x): no non-missing arguments to min; returning Inf

Warning in max(x): no non-missing arguments to max; returning -Inf

Warning: `stat_contour()`: Zero contours were generated

Warning in min(x): no non-missing arguments to min; returning Inf

Warning in max(x): no non-missing arguments to max; returning -Inf

Assemble figure 2

(g.mostlikely + coord_fixed(0.5))+
  plot_spacer() +
  ((g.fig2.acc / (g.fig2.corr + plot_layout(tag_level = 'new'))) + 
     plot_layout(heights = c(3, 1)) & 
     theme(plot.title = element_text(vjust = -8),
           legend.box.margin = margin(l=-5),
           strip.clip = "off",
           strip.text = element_text(size = 8),
           axis.text.x = element_text(angle = 45, hjust = 1))
  ) + 
  plot_layout(widths = c(1, 0.1, 2)) +
  plot_annotation(tag_levels = "a") & 
  theme(plot.tag.position = c(0, 1),
        plot.tag = element_text(vjust = 3)) -> g.fig2

g.fig2

if (save) {
  g.fig2 %>% ggsave(filename = "./figures/fig2_MattersArising.pdf", width = 8, height = 5.5, device = cairo_pdf)
}

Supplementary Figure 1

Prepare data

# Load subject-level data (used for clustering) and derive grand average (used for visualization) **without cropping**
tic()
df.supfig1.mean <- tibble()
df.supfig1.subj <- tibble()
for (manip in c("","_reord")) { 
  
  for (direction in c("rd_to_rd", "rd_to_mm", "rd_to_mp", "rd_to_or")) {
    tmp <- load_loop_subj(folder = "sounds",
                          condition = paste0(direction, manip),
                          crop = FALSE)
    
    df.supfig1.mean %<>% bind_rows(tmp[[2]] %>% mutate(manip = manip, direction = direction))
    df.supfig1.subj %<>% bind_rows(tibble(manip = manip, direction = direction, array = list(tmp[[1]])))
  }
}
toc() # ~1s

1.195 sec elapsed

Plot

# Extract the diagonal elements from the temporal generalization matrix
tmp <- filter(df.supfig1.subj, direction == "rd_to_rd", manip == "")[[1, "array"]][[1]]

tmp.diag <- matrix(nrow = dim(tmp)[3], ncol = dim(tmp)[1])
for (i in 1:dim(tmp)[3]) {
  tmp.diag[i,] <- tmp[,,i] %>% diag
}

# Plot
df.supfig1.mean %>% 
  filter(t_train == t_test,
         direction == "rd_to_rd",
         manip == "") %>% 
  mutate(sd = tmp.diag %>% apply(MARGIN = 2, FUN = sd) %>% `/`(sqrt(33))) %>% 
  filter(t_train >= -300, t_train <= 700) %>%
  
  ggplot(aes(x = t_train, y = accuracy)) +
  geom_hline(yintercept = 0.25, color = "grey50", linetype = 2) +
  geom_ribbon(aes(ymin = accuracy-sd, ymax = accuracy+sd), alpha = .3) +
  geom_line(size = 0.8) + 
  scale_x_continuous(breaks = scales::pretty_breaks(10), name = "Time (ms)") +
  scale_y_continuous(breaks = seq(.20, .40, .02), limits = c(.2,.4), name = "Accuracy") +
  coord_cartesian(expand = F)-> g.supfig1

g.supfig1

if (save) {
  g.supfig1 %>% ggsave(filename = "./figures/supfig1_MattersArising.pdf", width = 4, height = 2.5, device = cairo_pdf)
  g.supfig1 %>% ggsave(filename = "./figures/supfig1_MattersArising.png", width = 4, height = 2.5, device = ragg::agg_png(), dpi = 600)
}

Supplementary Figure 3

Functions & parameters

# Function to calculate and extract BF
getBF01 <- function(y, prior = "medium") {
  ttestBF(y, rscale = prior) %>% extractBF(onlybf = T)
}

# Color palette for BF
col.pal.bf <- rev(pals::brewer.piyg(7)); col.pal.bf[c(4)] <- "grey50";

Prepare data

tic()

## BF for EMPIRICAL null, DIFF
# --- extract array of subject-level correlation between accuracy and entropy 
X <- df.corr.subj[[3,"array"]][[1]]
# --- calculate BF for each cell of the TG matrix
bf.corr.diff.empir <- apply(X, c(1,2), getBF01, simplify = TRUE)


## BF for THEORETICAL null, DIFF
# --- extract array of subject-level correlation between accuracy and entropy 
X <- df.corr.maths.subj[[3,"array"]][[1]]
# --- calculate BF for each cell of the TG matrix
bf.corr.diff.maths <- apply(X, c(1,2), getBF01, simplify = TRUE)

## BF for MOST LIKELY null, DIFF
# --- extract array of subject-level correlation between accuracy and entropy 
X <- df.corr.ml.subj[[3,"array"]][[1]]
# --- calculate BF for each cell of the TG matrix
bf.corr.diff.ml <- apply(X, c(1,2), getBF01, simplify = TRUE)

toc() # ~1.5min

87.212 sec elapsed

## Post-processing
df.plot <- bind_rows(
  # --- aggregate data
  left_join(matrix_to_df2(bf.corr.diff.empir, "t_train", "t_test", "BF"), # BF results
            filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
    mutate(approach = "EMPIRICAL\nNULL"),
  left_join(matrix_to_df2(bf.corr.diff.maths, "t_train", "t_test", "BF"), # BF results
            filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
    mutate(approach = "THEORETICAL\nNULL"),
  left_join(matrix_to_df2(bf.corr.diff.ml, "t_train", "t_test", "BF"), # BF results
            filter(df.clusters.corr, manip == "original", stim == "sounds")) %>% # append correlation clusters
    mutate(approach = "MOST LIKELY")) %>% 
  
  # --- order factor levels 
  mutate(approach = factor(approach, levels = c("EMPIRICAL\nNULL", "THEORETICAL\nNULL", "MOST LIKELY"))) %>% 
           
  # --- convert indices to time
  mutate(t_test = (dt * (t_test - 1) - time_offset)) %>% 
  mutate(t_train = (dt * (t_train - 1))) %>% 
  # --- discretized version of BF
  mutate(BFd = cut(BF,
                   breaks = c(0, 1/30, 1/10, 1/3, 3, 10, 30, Inf),
                   labels = c("very strong in favor of null (BF < 1/30)",
                              "strong in favor of null (1/30 < BF < 1/10)",
                              "moderate in favor of null (1/10 < BF < 1/3)",
                              "inconclusive (1/3 < BF < 3)",
                              "moderate in favor of alt. (3 < BF < 10)",
                              "strong in favor of alt. (10 < BF < 30)",
                              "very strong in favor of alt. (BF > 30)")),
         .after = 3
         ) %>% 
  # --- transform continuous BF for easy symetric display on colorbar
  mutate(BF = case_when(BF < 1/3 ~ -1/BF,
                        between(BF, 1/3, 3) ~ NA,
                        BF > 3 ~ BF))

Joining with `by = join_by(t_train, t_test)`
Joining with `by = join_by(t_train, t_test)`
Joining with `by = join_by(t_train, t_test)`

df.plot %>% summary

    t_train        t_test           BF            
 Min.   :  0   Min.   :-670   Min.   :-5.000e+00  
 1st Qu.: 80   1st Qu.:-340   1st Qu.:-5.000e+00  
 Median :160   Median :   0   Median :-5.000e+00  
 Mean   :160   Mean   :   0   Mean   : 2.063e+06  
 3rd Qu.:240   3rd Qu.: 340   3rd Qu.:-4.000e+00  
 Max.   :320   Max.   : 670   Max.   : 1.539e+10  
                              NA's   :4029        
                                          BFd           signif      
 very strong in favor of null (BF < 1/30)   :   0   Min.   :0.0000  
 strong in favor of null (1/30 < BF < 1/10) :   0   1st Qu.:0.0000  
 moderate in favor of null (1/10 < BF < 1/3):7945   Median :0.0000  
 inconclusive (1/3 < BF < 3)                :4029   Mean   :0.2496  
 moderate in favor of alt. (3 < BF < 10)    : 569   3rd Qu.:0.0000  
 strong in favor of alt. (10 < BF < 30)     : 264   Max.   :1.0000  
 very strong in favor of alt. (BF > 30)     : 558                   
    manip               stim                      approach     
 Length:13365       Length:13365       EMPIRICAL\nNULL  :4455  
 Class :character   Class :character   THEORETICAL\nNULL:4455  
 Mode  :character   Mode  :character   MOST LIKELY      :4455  
                                                               
                                                               
                                                               
                                                               

Plot

lines.v.pos <- 330*seq(round(min(df.plot$t_test)/333), round(max(df.plot$t_test)/333), 1)
lines.v.pos <- lines.v.pos[lines.v.pos!=0]
lines.h.pos <- 330*seq(round(min(df.plot$t_train)/333), round(max(df.plot$t_train)/333), 1)
lines.h.pos <- lines.h.pos[lines.h.pos!=0]


df.plot %>%
  ggplot(aes(x = t_test, y = t_train)) +
  facet_grid(approach ~ .) +
  geom_tile(aes(fill = BFd), show.legend = TRUE) +
  
  geom_vline(xintercept = lines.v.pos, linetype = 2, color = "black") +
  geom_hline(yintercept = lines.h.pos, linetype = 2, color = "black") +
  geom_vline(xintercept = 0, color = "black") +
  geom_hline(yintercept = 0, color = "black") +
  
  geom_contour(aes(z = signif), size = 0.2, color = col.clusters) +
  
  scale_x_continuous(breaks = scales::pretty_breaks(12),
                     expand = c(0,0)) +
  scale_y_continuous(breaks = scales::pretty_breaks(4),
                     expand = c(0,0)) +
  scale_fill_manual(values = rev(col.pal.bf),
                    labels = levels(df.plot$BFd),
                    drop = FALSE) +
  
  labs(fill = expression(BF[10])) +
  coord_equal() +
  
  labs(x = "Test time (ms)", y = "Train time (ms)") -> g.supfig3

g.supfig3

if (save) {
  g.supfig3 %>% ggsave(filename = "./figures/supfig3_MattersArising.pdf", width = 10, height = 3, device = cairo_pdf)
  g.supfig3 %>% ggsave(filename = "./figures/supfig3_MattersArising.png", width = 10, height = 3, device = ragg::agg_png(), dpi = 600)
}

Supplementary Figure 4

mat.C.all <- array(dim = c(34,34,4,4,n.subj))
i <- 1

for (s in list.subj) {
    
    # --- load the confusion matrix
    tmp.file <- os$path$join(path.root, s, "rd_to_rd_confmats.npz") %>%
      np$load()
    
    # --- crop around S0 (decoded stimulus)
    mat.C.all[,,,,i] <- tmp.file$f[["arr_0"]][,71:104,71:104,,] %>%
      apply(c(2,3,4,5), mean)
    
    i <- i+1
}

mat.C.mean <- mat.C.all %>% apply(c(1,2,3,4), mean)

    
df.confusion <- matrix_to_df2(mat.C.mean[22,22,,], "from", "to", "p") %>% 
  mutate(label = paste0("a[",from,to,"]"))
df.confusion %>% 
  ggplot(aes(x = to, y = from, fill = p)) +
  geom_tile(color = "firebrick3") + 
  geom_text(aes(label = label), parse = T, color = "firebrick3") +
  scale_x_continuous(name = "predicted class", position = "top") +
  scale_y_reverse(name = "true class") +
  scale_fill_stepsn(name = "p",
                    colors = pals::brewer.greys(100),
                      # low = "white", high = "black",
                      labels = ~str_replace(., "^0.","."),
                      breaks = scales::pretty_breaks(6),
                      limits = c(.20,.30), oob = scales::squish) +
  coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
              expand = F, clip = "off") -> g.matC

matrix_to_df2(mat.T$MP, "from", "to", "p") %>% 
  mutate(label = paste0("a[",from,to,"]")) %>% 
  ggplot(aes(x = to, y = from, fill = p)) +
  geom_tile(color = "firebrick3") + 
  geom_text(aes(label = label), parse = T, color = "firebrick3") +
  scale_x_continuous(name = "to", position = "top") +
  scale_y_reverse(name = "from") +
  scale_fill_stepsn(name = "p",
                    colors = pals::brewer.greys(100),
                      labels = ~str_replace(., "^0.","."),
                      breaks = scales::pretty_breaks(6),
                      limits = c(.0,.50), oob = scales::squish) +
  coord_equal(xlim = c(0.5,4.5), ylim = rev(c(0.5,4.5)),
              expand = F, clip = "off") -> g.matT

 
((g.matT + labs(caption = "Transition matrix\n(midplus (MP) sequence)") +
    guides(fill = guide_colorbar(barwidth = 0.5, ticks = F))) +
(g.matC + labs(caption = "Confusion matrix\n(group average, t=210ms)") +
    guides(fill = guide_colorbar(barwidth = 0.5, ticks = F))) &
  theme(plot.caption = element_text(size = 10),
        axis.title = element_text(face = "italic"),
        legend.box.spacing = unit(2,"mm"),
        )
  ) +
  plot_annotation(tag_levels = "a") -> g.supfig4

g.supfig4

if (save) {
  g.supfig4 %>% ggsave(filename = "./figures/supfig4_MattersArising.pdf", width = 7, height = 2.5, device = cairo_pdf)
  g.supfig4 %>% ggsave(filename = "./figures/supfig4_MattersArising.png", width = 7, height = 2.5, device = ragg::agg_png(), dpi = 600)
}

References

Demarchi, Gianpaolo, Gaëtan Sanchez, and Nathan Weisz. 2019. “Automatic and Feature-Specific Prediction-Related Neural Activity in the Human Auditory System.” Nature Communications 10 (1): 1–11. https://doi.org/10.1038/s41467-019-11440-1.